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Consider a dynamical system f : M. ^ AA, on a manifold A4 endowed with the 
Lebesgue measure. Assume that there is a unique /-ergodic probability measure p so 
that limT^oo Y.I=i 4>{f^ i^)) = I 'pdp for continuous functions and for Lebesgue- 
almost every x & M.. Then, p is called the SRB measure or physical measure (see 
[21] for a discussion on these two terminologies). Assume that when / is perturbed 
in its neighborhood, the average value p{A) of a given smooth observable A varies 
differentiably with /. Then the derivative operator is called the linear response function. 
In the case of uniformly hyperbolic systems, and provided that / and its perturbations 
are sufficiently smooth, then the linear response is given by [T9] : 

oo 

5p{A) = Y.p{X.V{Aon), (1) 

?i=0 

where X is a (time- independent) perturbation. Ruelle has indeed shown that this series 
is convergent for suitable perturbations. In the case of a time-dependent perturbation 
Xt this formula extends formally to: 

oo 

5pt{A) = Y.{K^Xt_^)A. (2) 

The operator Hn : Ti. ^ B* , where 7i is a suitable space of vector fields and B* the dual 
of a suitable space of functions, is defined by: 

(/€„Xi_„)A = p(Xi_„.V(Ao/")); n>0, (3) 

while, for n < 0, (/t„Xt_„)y4 = (causality). It is called the "response function". In this 
case, one recovers the standard definition of the linear response in non equilibrium 
statistical mechanics (([2]) is a convolution product), and classical results such as 
fiuctuation-dissipation relations or Onsager theory can be recovered [I9]. Moreover, 
this result does not require the assumption of closeness to equilibrium. 

However, one can find examples of non uniformly hyperbolic systems where the 
series ([1]) does not converge, even for maps of the interval with absolutely continuous 
invariant measures [18]. In this case, one considers the complex susceptibility function, 
x(a;) = {ki^X)A where k^^ is the Fourier transform of the response function. It is also 
convenient to use the variable A = e*"^. Then, x{^) = ^('^)- Under some conditions on 
X and / (interval maps), Ruelle proved that \1/(A) extends to a meromorphic extension 
with no pole at A = 1. He also showed that this function has a pole inside the 
unit disk |A| < 1 (i.e. in the upper-half complex plane for the frequencies u) and 
other poles, corresponding to the eigenvalues of the Perron-Frobenius operator (Ruelle- 
Pollicott resonances [121 E])) outside the unit disk (in the lower-half complex plane). 
The existence of a pole of the complex susceptibility in the upper-half complex plane 
must not be misinterpreted. This is not a "violation of causality" , which has no meaning 
in a dynamical system, causal by construction. This simply expresses that an arbitrary 
small perturbation of / does not give a response proportional to the perturbation. 
According to Ruelle: "this might be expressed by saying that p is not linearly stable" . 
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This result is quite interesting and intriguing for a physicist and two natural 
questions arise. Does this property exist in larger dimensional chaotic dynamical 
systems ? And is it possible to have an experimental/numerical characterization of 
the effects of this pole ? The first question is difficult to address from a mathematical 
point of view and is still unresolved even for classical examples such as the Henon 
map. The second one is the main concern of the present paper. We indeed propose a 
numerical experimentation protocol allowing to numerically approximate the complex 
susceptibility and to investigate the numerically observable effects of having a pole 
in the upper-half complex plane. 

The paper is organized as follows. After a brief summary of Ruelle's results for the 
logistic map (section [1]), we describe our numerical methods in section [2] and apply it to 
the logistic map where one can use the guidelines of Ruelle's mathematical results as a 
validation (section [3]). Then, we apply it to the Henon map for the standard parameter 
values and conjecture that this map could also exhibit a pole in the upper-half complex 
plane (section H]). 



1. Reminder of Ruelle's results. 

The dynamical system defined by the logistic map: 

/(x) = l-2x^ XG[-1,1], (4) 
has an absolutely continuous ergodic measure with density : 



The linear response ([I]) reads, in this case: 

d 

n=0' 



SpiA) = E / X{x)-A{rx)p{dx), (6) 
while the complex susceptibility reads: 

oo „ 

xicu) = M/(A) = E ^"7 A'{rx){rnx)x{x)p{dx) (?) 



n=0 



with A = e*"^. Ruelle has shown the following results Consider the variable change 
X = sin(|?/) = zu{y) which maps / on g{y) = 1 — 2\y\ [tent map]. Then the density p is 
mapped onto the density cro{y) = |. In this variable, we have : 

vl/(A) = EA"/ iC-Y){s)B'is)ds (8) 

n=0 -^^^ 

where B = Aou7, Y{y) = My)^^^ and Co^y) = $ (^) - $ (i^). 

Consider first a square integrable holomorphic perturbation X = Xq that vanishes 
at {—1,1}. Then, the series ([8]) has a meromorphic extension in C and has no pole 
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inside the unit disc. Moreover, it is holomorphic at A = 1 (resp. u = ). Indeed, one 
has, using integration by parts: 

^(A) = -EA"/ {C-Y^){s)B{s)ds, (9) 

n=0 

where C is the Perron-Frobenius operator: 

mv) = i* f^l + i* C-^) (10) 



2 V 2 / 2 V 2 

In the space of analytic functions, and in the basis of monomials {y^}, m = 
0, 1, 2, . . . this operator is represented by an infinite triangular matrix. The eigenvalues 
are on the diagonal. They are of the form 1/4", n = 0,1.... The corresponding 
right eigenf unctions are (Bernoulli) polynomials such that cTo{y) = \ (invariant density 
corresponding to the eigenvalue 1); 

aM = -\-'^y + y^ (11) 

(eigenvalue 1), etc... This simple argument can be found in [TU] pp 123-125. This result 
is also a direct consequence of [3]. 

Therefore, for a perturbation Yq{ii) = Y^'^=iCkcrk{y), and for an observable A such 
that the integral j\YQ(s)B{s)ds is finite, one finds easily that the susceptibility has 
poles \k^,k > 1, outside the unit disk. For example, if Vo(^) ~ '^liv): the series ([9]) 
reads: 



n=l 



^W = -E|4) I ^^iis)B{s)ds (12) 



-1 



This case is used as a benchmark for the numerical procedure developed in the next 
section. The corresponding perturbation Xq reads, in the variable x: 

2 . _ 4 



Xo{x) = ttVT^^ 
where Xo(±l) = 0. 



arcsin(a;) arcsin^(a:;) H arcsin^(x) 

3tt n-^ Stt'^ 



(13) 



Ruelle has also shown that there is a meromorphic function : 

^_(z) = -^ + w{z), (14) 



where w{z) is of order 0{z + 1), and where $_(±1) = such that Co^~ = y f'i~^)^- 
(= 2$_ in our case). The function w in (flU) is therefore a solution of: 



w is not an analytic function. It can be obtained via the recursion: 

\ -w(l-2z)-^; 2>0; 
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reminiscent of the Tagaki function discussed by Gaspard in the context of the multi- 
baker map [To]. Returning to the variable x, the function $„ corresponds to : 



X(x) = nVT^ 



x^ 



1 ,2 

-2 . , . . ^ + w[ — arcsin(x) + 1) 



■ arcsin(x) + 1 tt 



:i7) 



The corresponding series ^_(A) = E^=o ^" /p(a;)^'(/"a;)(/")'(x)X_(x)rfx, 
converges if |A|w|/'(— 1)| < 1, and, in its domain of convergence, has the form 



1 ^1^1/' ( 1)1 ^^^^^ ^ constant depending on A. The complex susceptibility has 

therefore a pole at A = ^^^^ = |, inside the unit disk (resp. u = 'jlog(2)). Therefore, 
the series diverges for real frequencies uj. 

2. The numerical method. 

The numerical method used here has been introduced by Reick in pfS] and independently 
in [71 El [9]. It is purely heuristic and there is no mathematical proof of its validity, though 
a discussion with some analytic developments can be found in [13] and [9] . Thus, in this 
section, we shall use formal expressions. 

Consider a time-dependent perturbation of / of typt 

xt+i = f[xt + eXtixt)] . (18) 
The variation of the average value of A, at time t, is formally given by: 

ApM) = e E / A'{f--x){f-n'{x)X4x)p{dx) + 0{e'). (19) 

n=— oo 

which corresponds to ([2]) to the linear order and up to a factor e. Note that the term 
O(e^) is not uniform in t. 

Consider now perturbations of type Xt{x) = X(a;)e~*'^* where u is real. Then: 

Apt{A)=ee-'-'x{u^) + 0{e'), (20) 

where x{^)y the complex susceptibility for X{x), is given by eq. ([7]). Then, at the linear 
order in e: 

e 

Since xi^) does not depend on time one may write: 

x{^) = ^j:^Pt{A)e-' (22) 
^-^ t=i 

I Note that this way of perturbing avoids the introduction of /^^ in the computation of the response 
function as found e.g. in [19j . It allows one to write convolution formulae and the Fourier transform 
in the standard form they have for continuous time systems. This was the main motivation for this 
choice in ^ and we pursue along these lines in the present paper. 
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The time-dependent SRB state pt involves an average over initial conditions in the 
distant past [19]. One can argue (see [13] Appendix A) that the above average over t 
makes the average over initial conditions unnecessary. The idea is therefore to replace 
/S.pt{A) by A{x' it)) — A{x{t)) where x'{t) is a typical trajectory of the perturbed system, 
with a perturbation X(x)e"*'^* and x{t) is a typical trajectory of the unperturbed system. 
Then, one computes: 

^t(^, x) = ^ E [A{x\t)) - A{xm e-* (23) 
t=i 

The hope is that S't(co',x) x{^) as T oo, where the limit is independent of 
x(0). There is no mathematical proof of the above statement. Note that T must be 
taken so large that roundoff errors play their role of selecting the limit to be the SRB 
state. On practical grounds, as discussed in [T3], this procedure is numerically rehable 
provided that eT is sufficiently large and uT » 1 (basically one uses eTu >> 1). These 
conditions have been checked in the simulations performed in the present paper. Since 
St{uj,x) is expected to give a fairly good approximation of x{^)^ the Fourier transform 
of the linear response, when T is sufficiently large, one also expects that the inverse 
Fourier series of St{uj,x), called R{t) (where we drop the T and x for simplicity) gives 
a fairly good approximation of the linear response ([2]). 

Remarks 

(i) The quantity obtained by this procedure is, stricto sensu, not the linear response 
but the total response since one can apply it for arbitrary large e. As a test, one 
must a posteriori check that the susceptibility is independent of e on a certain range 
of small e values [e.g. it does not vary if one replaces e by 2e]. 

(ii) In the case of the logistic map, any typical trajectory of approaches the points 
±1 within a distance of order e with a characteristic time of order In a numerical 
simulation where e is small but finite this arises often, especially if one respects the 
condition eTuj >> 1. However, in this perturbation X{x)e can push 
the trajectory out of the interval [—1, 1] leading to an exponential divergence. To 
avoid this effect we have used a perturbation X{x,t) = X{x){l + e^*"^*) where X 
is positive about x = — 1 and negative about x = +1, so that the perturbation is 
always directed inside the interval [—1, 1] whenever x = ±1. 

3. Application to the logistic map. 

3.1. A benchmark. 

Let us first consider the case where the complex susceptibility is well defined on the real 
axis (for the frequency lu) by the series (JTj). More precisely, we consider a perturbation, 
corresponding to the Ruelle-PoUicott resonance |. In the variable x this perturbation 
is given by eq. (fT3l) . For the observable A one can use any choice such that the integral 
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j\ ai{s)B{s)ds in ([I2]) is finite. We choose A = Xq. The complex susceptibihty reads 
^(A) = where C > is a constant. 

In the figure [1^ the complex susceptibility, computed with this algorithm and 
compared to the theoretical value, is drawn. To generate the real and imaginary 
parts of the susceptibility we run two independent simulations. Starting from the same 
mother trajectory, we apply in the first case a perturbation ecos(u;t) and compute 
the corresponding susceptibility, which gives Re{x{^))- In the second case we apply 
a perturbation — esin(ti;t), which gives Im{x{^))- Note that the perturbations are 
therefore real and the trajectories stay on the real axis. The amplitude of the 
perturbation was fixed to e = 10~^, 5 x 10^'^, 10~^. Averages are performed by computing 
the time average St{u},x) (l23|) for initial conditions Xn,n = 1...N randomly 
chosen and then computing Y.n=i Sri^yXn)- This procedure allows us to reduce the 
numerical noise and to compute error bars. We took T = 10^, = 10 for e = 10~^; 
T = 2.10^A^ = 20 for e = 5 X lO^^ and T = 10^ A^ = 10 for e = 10^^. Note that 
numerical noise is large and that the resonance curve is very fiat, requiring to have 
small error bars. This requires an average over a very long time T, roughly given by the 
condition eTu >> 1, and limits the range of e values that one can reach. Note that the 
condition eTuj >> 1 is always violated as u approaches zero, whatever the (finite) value 
of T. This explains the discrepancy observed for small frequencies. If one excludes this, 
the agreement is quite good. Note also that the experimental curve has a magnitude 
that does not depend of e, in this range of e values, as expected. 

The (approximate) linear response R{t) can be easily obtained by computing the 
inverse Fourier series of x{^) (s-g- by a direct summation). It is drawn in Figure [lb 
(in log scale for the y axis). Note that up to numerical noise R{t) is real, as expected. 
One observes an exponential decay very close to the theoretically expected decay 4~*, 
corresponding to the mixing rate given by the first Ruelle-PoUicott resonance (and also 
to the characteristic decay rate toward equilibrium, in agreement with the fiuctuation- 
dissipation theorem). Note that the range of validity for the interpolation is very thin 
(~ [0,4]). Indeed, rapidly the perturbation becomes so weak that one measures only 
the numerical noise. 

3. 2. Perturbation X_ . 

As we saw, the perturbation (JTTl) leads to a diverging series (171) for real frequencies. 
However, when computing the expression fl23l) one is not dealing with a series, but with 
a finite sum, which diverges as T grows. The divergence rate provides useful hints on 
the presence of the pole. 

Consider indeed the formal expression (fT9l) of Apt{A) in the case of the logistic map. 
jt-n Y^^^ t — n zeros and oscillates rapidly (for large t — n) with a period ~ 4"~*. One 
can decompose the integral p{x)A' {f^'^x){f^^'^y {x)X{x)dx over intervals delimited 
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Figure 1. (a) Susceptibility for a perturbation and an observable also given by 
(fT3|l . for e = 10~^,5 x 10~^,10~'^. (b) Corresponding approximate linear response 
[denoted by R(t)], in log scale and theoretically expected curve y = a — log{4)t. 



by the zeros of The density p has strong variations about ±1 and small variations 
in the central part. The contribution of the first interval (containing —1) is given, for 
large t — n, by : 



r^^""* p{x)A'{r-x){f-)'{x)x{x)dx 

where x = —1 + u. If 74'(/*""(— 1 + u))X{—l + u) does not vanish about m = (x = —1) 
and has small variations on this interval, this contribution is ~ 4*"*^ in -i=du = 2*""^, 
which diverge^j when t—n oo, with a rate exactly given by the pole A = 2- Therefore, 
another good benchmark of our method, in the case where the complex susceptibility 
has a pole in the upper-half plane, consists in computing fl23|) . and taking the inverse 
Fourier series, providing in this way an approximation for Apt(^)- One should see then 
the exponential divergence with a rate log(2) (which is exactly the value of the Lyapunov 
exponent A). 

However, there is another effect that must be taken into account in the numerics. 
When the amplitude of the response, ~ e2*, is small one computes numerically the 

§ If one considers now the contribution of the "bulk" [inner intervals] to the integral ([24|l one can 
figure out that it does not diverge, essentially because the derivative (/*^")'(2;) has alternating signs 
and because the density p has small variations in the bulk. More precisely, the contribution of the bulk 
is provided by a decomposition on the Ruelle-PoUicott eigenfunctions (7„ and decays exponentially due 
to mixing. Thus, the existence of a pole in the upper-half complex plane is due to a boundary effect 
obtained when the perturbation weights the points {—1, 1} (where the density diverges). 
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first term of tlie Taylor expansion (1191) . But when e2* becomes too large the numerical 
computation includes also the nonlinear terms and one computes in fact the complete 
susceptibility, including nonlinear terms. This one does not diverge because the trajec- 
tories of the perturbed system remain in the interval [—1,1], but it does not depend 
linearly on e. 

Therefore, one expects the following. There is a time cut-off tm(e) ~ ~Iog(^ 
beyond which one does not compute the linear response but the complete response, 
including nonlinear effects which saturate the growth of the perturbation. Then one 
should observe an exponential growth 2* up to tm(e) in the linear response of the 
form 2*i7[o,t,„(e)](^); where Hj{) is the characteristic function of the interval /. The 
corresponding susceptibility is: 

which reads, using the frequency u: 

(l~2cos(cj))(l-2*'" cos(a;f,n))+2*'"+^ sin(Lj) sinji^tm) 
5-4cos(aj) 

, • 2sin(a;)(l-2''" cos{ujtm))-2*"^ sm(a;f^)(l-2 cos(lj))) K^^^J 
5-4cos(lj) 

Therefore, it exhibits oscillations due to the cut off tm- For longer time, mixing should 
lead to an exponential decay of the response. 

To check this we have first computed (123!) for the perturbation ( [T7|) and for 
the observable A{x) = x — ^ (the derivative of A gives a maximal value for the 
term A'{f^~"'{—1 + u)) in equation fl2^ ). The function X_ has been computed with 
the recursion up to the order 4. The e values 10"^ 10-^10-^(T = 10^, N = 
400),10-5(r = 10^, N = m),10-%T = 4 x 10^, N = leOO) have been considered (but 
only 10~^, 10~^, 10~^ are represented in Fig. [2K,b for the legibility of the figure.). We 
have used a fit procedure to compare the experimental data with the form fl25l) . We 
have also computed the approximate linear response. The results are presented in Fig. 
[3i,b. 

One observes indeed a clear dependency with e. Moreover, the real and imaginary 
part exhibit the expected oscillations due to the cut-off, with a perfect agreement with 
the form fl25l) . Taking the inverse Fourier series one obtains the linear response in Fig. 
[2]d. One sees clearly the exponential growth 2* and the e dependent cut-off. Performing 
a fit in log scale one obtains (fig. [3]) an exponential increase with a rate 0.65 very close 
to the expect value log{2) = 0.693. After this there is an exponential decay with an 
approximate rate —0.24. (We have only represented the fit for e = 10~^ but one sees 
easily in Figl3] that the decay rate is similar for e = 10~^,10~^). This rate is slower 
than the first Ruelle Pollicott resonance — log(4) = —1.386. It might be that we are 
observing a crossover regime where exponential amplification and exponential mixing 
are competing. 
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Figure 2. (a) Susceptibility for tlie perturbation where w was computed up 

to order 4, and observable A{x) — x~ The parameter e takes the value 10^^, 10^^. 
In full line are drawn the fitting curves obtained from (|25p . (b)Linear response. 




Figure 3. Fit of the linear response corresponding to Fig. [2}d. 



4. The Henon map. 

The Henon map: 

= 1 - axj + yt 




with a = 1.4, b = 0.3 has a strange attractor with a fractal structure. Moreover, 
numerical estimates of the largest Lyapunov exponent, using e.g. Eckmann-Ruelle 
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algorithm [S] gives a value ~ 0.422. But the Henon map is not uniformly hyperbolic. 
There are points with tangencies between stable and unstable manifolds. These points 
may be responsible for phenomena analogous to those described in the previous section, 
for specific perturbations weighting those points. From this observation, it is conjectured 
that there may also exist a pole in the upper-half complex plane (and possibly more 
complex singularities |20]). However, there is no mathematical result for this and the 
form of the perturbations/observable leading to such a singularity is not known. 

Using the same method as in the previous section, a natural empirical approach 
consists in applying a perturbation ee*"^* in one of the directions (x, y) and investigating 
the effect on the linear response for one of the variables {x,y). As an example, 
we have numerically computed the response of the observable A{x,y) = x to a 
perturbation ee*"^* in the direction y. Denote by Xxy{^) and Rxyif) the corresponding 
susceptibility and approximate response. They are drawn in Fig. |l^,b respectively, 
for e = 10-3(T = 10^ iV = 100),10-^(T = 10^ iV = 100), 2 x 10-^(T = 10^ = 
400), 10~^(T = 10^, = 1600). Note that essentially the same curves are obtained, up 
to a phase factor, when perturbing direction y and computing the response of y. In 
particular \Xxy{'^) \ = 1x^(^)1- 




■3 -2,5 -2 -1,5 -1 -0,5 0,5 1 1,5 2 2,5 3 
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Figure 4. (a) Susceptibility for a constant perturbation X = e in the direction y for 
the observable A{x,y) = x. (Jo) Corresponding linear response. 



One observes similar effects as in the previous section. The amplitude of the 
complex susceptibility increases with e (and this effect is not stabilized, neither by 
increasing the time T of average well beyond the criterion eooT >> 1, nor by increasing 
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the number of sample trajectories). The hnear response has three parts. For short 
times [t < 20) there is a small bump for the modulus of the response: its height does 
not seem to depend on e. The perturbation/observables have a nonzero projection onto 
the Ruelle-Pollicott modes. It might be that this bump comes from this part, having 
a well behaved linear response. The second part exhibits an exponential increase with 
a time cut-off depending on e. This is very similar to the observations made for the 
logistic map: exponential amplification of the linear response until nonlinearities induce 
saturation of this effect. Then, mixing leads to an exponential decay (third part of the 
curve). 

The exponential instability is well fitted by the curve e°e''^* with a = —6.80 ± 0.87 
and A = 0.455 ± 0.038 (see Fig. [5]). Note that the positive Lyapunov exponent 
of the Henon map is 0.42(2) for these values of a, b. The decay in the third part 
is well interpolated by a sum of exponential terms. The dominant term is /3e'''* 
with 7 = —0.129 ± 0.015. The value of 7 is in agreement with the exponential 
decay of the correlation function Cxyif) =< f^{x)y > — < x >< y >, where 
< g >= J g{X)p{dX), X = {x,y) and p is the SRB measure ( Fig. [5]). Thus the 
last part apparently obeys the fluctuation-dissipation theorem. Correlation function 
has been computed by the standard Fast Fourier Transform method (see [11], chapter 
13.2). 

Let us also remark that there are no thin peak in the susceptibility (the frequency 
resolution is 0.00612) and the thickness of the peaks does not change with e. Since the 
widtlJJJ] of the resonance peak is given by the imaginary part of the pole, this suggests 
that the imaginary part of the poles is bounded away from zero. This leads us to 
conjecture that there is no pole on the real axis (for the frequency 00). This conjecture 
is compatibl^ with the available results (in particular exponential decay of correlations 
for many parameters a and b (see [H O [231, [22])). 

5. Conclusion. 

There is currently an intensive research activity in mathematics dealing with linear 
response theory for interval maps (see e.g. [H [2] and reference therein). In the case 
of non uniformly hyperbolic maps one can find examples where the susceptibility has 
an arbitrary large number of poles in the upper-half complex plane [21] or where the 
average of a smooth function with respect to the SRB measure of the perturbed system 
/ + eX is not even Lipschitz at e = [1] . These results are intriguing for a physicist and 
raise some natural questions: can we measure some effects related to these situations? 
What is the "genericity" of these examples? Do these effects arise in dynamical systems 

II Consider the complex susceptibility in the neighbourhood of a simple pole ujq — xq + iyQ. Then 
x(w) ~ ui-ujo • A_ssume that this pole is close to the real axis. Set ui = x + iy. On the real axis 
|x(w)p ^ (^x~xl)^+y'^ • '^^^ modulus is maximal at x = a;o: ^^id has value ^ (resonance peak). The 
width of the resonance peak is 2\xi — xo\ where xi is such that [x(xi)| = Thus, the width is 2yo- 
% I thank one of the reviewers for this remark. 
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Figure 5. (a) Fit of tlic linear response \Rxy{t)\, e = 10 ^ for the Henon map. (b) 
Fit of the correlation function C'xy(t) for the Henon map. 

modeling physical systems ? do they occur in larger dimensions? Answering these 
questions at the level of rigor of mathematics is probably a too formidable task at the 
current state of research in this field. One hope is that some hints can be provided by 
numerical approaches. 

In this spirit, we have presented in this work a numerical procedure allowing to 
compute the complex susceptibility/linear response. For this, we use an approximation 
of the complex susceptibility ([7]) by a finite sum, which can be easily computed. In the 
cases where the series ([7]) converges this sum gives a fairly good approximation. When 
there is a pole in the upper-half plane a correct treatment of this sum allows us to 
recover the pole location. Applying this procedure to the Henon model we conjecture 
that such a pole could also be present. 

There are several possible developments of the present work. On mathematical 
grounds it could be interesting to apply this procedure to the cases discussed in e.g. 
pTl [H, [2] that can be addressed by rigorous methods. (In particular the cases in [1], [2] 
should provide examples where the linear response does not decay exponentially, but 
does not grow exponentially either.) An intermediate case between these mathematically 
tractable cases and the difficult Henon model would be CoUet-Eckmann maps. Another 
question is: "what is there beyond the pole" ? In some sense, the pole in the upper- 
half complex plane may "hide" more complex singularities lying behind it such as cuts. 
Is it possible to have any (numerical) idea of which type of singularities are there ? 
A natural way of doing this is to remove the exponential instability by multiplying 
the perturbation by a damping factor (or equivalently to use a complex frequency uj). 
This is natural from a mathematical point of view, but tricky in the numerics, because 
the damping factor is either smaller than the exponential instability (and the response 
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grows rapidly inducing nonlinear effects, as we saw) or it is bigger (and the perturbation 
becomes rapidly numerically 0, then we are measuring short time transients). One has 
thus to make small variations of the damping factor around the pole and look at the 
changes in the susceptibility/response curve. 

On more physical grounds, a natural continuation of the present work, could be to 
consider a lattice of coupled logistic maps, to apply an harmonic perturbation at some 
point and look at the induced effects, in the spirit of this work. This would be one step 
towards the investigation the effects of large dimension ( "thermodynamic limit" ) on the 
singularity induced by the microscopic dynamics. 
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